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Abstract. Quasiperiodic oscillations and shape-transformations of higher-order 
bright solitons in nonlinear nonlocal media have been frequently observed in recent 
years, however, the origin of these phenomena was never completely elucidated. In this 
paper, we perform a linear stability analysis of these higher-order solitons by solving 
the Bogoliubov-de Gennes equations. This enables us to understand the emergence of 
a new oscillatory state as a growing unstable mode of a higher-order soliton. Using 
dynamically important states as a basis, we provide low-dimensional visualizations of 
the dynamics and identify quasiperiodic and homoclinic orbits, linking the latter to 
shape-transformations. 
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1. Introduction 

Bright solitons are particle-like nonlinear localized waves , that keep their form while 
evolving due to a compensation of diffraction or dispersion of the medium by the 
nonlinear self- induced modification of the medium [TJ. Usually, solitons are studied in 
systems exhibiting local nonlinearities, where the guiding properties of the medium at a 
particular point in space depend solely on the wave intensity at that particular point [2] . 
Here, we consider nonlocal nonlinearities, i.e. situations in which the nonlinear response 
of the medium at a point depends on the wave intensity in a certain neighborhood of 
that point, where the extent of this neighborhood is referred to as degree of nonlocality. 
Nonlocal nonlinearities are ubiquitous in nature, for example, when the nonlinearity is 
associated with some sort of transport process, such as heat conduction in media with 
thermal response [3H5], diffusion of charge carriers [6j[T] or atoms/molecules in atomic 
vapors [8,9j. Nonlinearities are also nonlocal in case of long-range interaction of atoms in 
Bose-Einstein condensates (BEC), such as in case of dipolar BEC [T0HT3] or BEC with 
Rydberg-mediated [T^[T5] interactions. In addition, long-range interactions of molecules 
in nematic liquid crystals also result in nonlocal nonlinearities [T6HT9] . 

The balance between diffraction and nonlinearity may lead to stable solitons 
withstanding even strong perturbations. In particular, it has been shown, that nonlocal 
nonlinearities crucially modify stability properties of localized waves. With respect to 
bright solitons, they lead to a much more robust evolution as compared to its local 
counterpart (20j[2TJ. This is due to the fact, that nonlocality acts like a filter by 
averaging or smoothing-out effect on perturbations which would otherwise grow in case 
of local response of the medium [22]. For example, higher- dimensional solitons would 
collapse for systems exhibiting local nonlinearities, whereas they can be stabilized by 
nonlocality [23H25]. 

In this work, we investigate the linear stability and nonlinear dynamics of higher- 
order solitons. In particular, we study the quadrupole soliton Q and the second-order 
radial soliton R 2 (a hump with a ring), as sketched in Fig. [TJ For those solutions, 
a quasiperiodic shape transformation between states of different symmetries has been 
observed recently in [26l[27J. However, a complete understanding of this spectacular 
phenomenon is still missing. One difficulty in the analysis of the shape transformations 
is that they cannot be described solely in terms of linear perturbation because they are 
not small (27]. Nevertheless, here we show that in spite of the fact that we are dealing 
with a highly nonlinear phenomenon, deeper understanding can be gained from the 
linear stability analysis of the corresponding Bogoliubov-de Gennes (BdG) equations. 
In other words, solutions of the linear stability analysis of the solitons are used to 
describe wave dynamics in the neighborhood of a soliton solution. Moreover, in order to 
fully understand nonlinear dynamics, we employ and further develop techniques recently 
introduced in dynamical systems studies of dissipative partial differential equations 
(PDE) [28l[29]. These methods employ projection of PDE solutions from a functional 
infinite space onto a finite number of important physical states or dynamically relevant 
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Figure 1. Two particular soliton solutions: a) quadrupole soliton Q and b) second- 
order radial soliton i?2- Both soliton profiles can be chosen real without loss of 
generality. The lower plane shows the modulus square depicted in color scale of the 
two solitons. 



directions. Here, the relevant directions are mainly the unstable and stable internal 
modes of the solitons. The introduction of these low- dimensional projections will 
allow us to interpret the non-periodic soliton oscillations as indication of homoclinic 
connections. Moreover, we are able to understand how different solutions, including 
quasiperiodic oscillations, are organized by this homoclinic connection. The same 
analysis should also work for a larger variety of higher-order solitons of this nonlocal 
system, such as those presented e.g. in [26|[2T] . 

The paper is organized as follows: In Sec. [2} we introduce the governing equations 
of motion. In Sec. [31 we solve the BdG equation to find the internal modes of 
the quadrupole soliton Q as well as the second-order radial soliton R 2 . In Sec. IU 
we discuss nonlinear soliton propagation, introduce low- dimensional projections and 
study homoclinic and quasiperiodic trajectories in this representation. Finally, we will 
conclude in Sec. El 

2. Model equations 

The underlying model equation for our subsequent considerations is the nonlocal 
nonlinear Schodinger equation (NLS) 



where A = d xx + d yy denotes the transverse Laplacian. Depending on the actual 
context, \ip(r, t)\ 2 can be identified with either the intensity of an optical beam in scalar, 
paraxial approximation, or the density of a two-dimensional BEC within mean field 
approximation. The nonlinearity 9 is given by the convolution integral 



where the kernel K is determined by the physical system under investigation, and 
r = (x,y). If K(r) = K(\r\), then Eq. (CQ) is invariant under rotation and the angular 



id t ip + Aip + dip = 0, 



(1) 




(2) 
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momentum is conserved. This is the case here, as we consider the Gaussian nonlocal 
model, for which quasiperiodic oscillations have been originally observed [26] [27]: 

K(r) = e" r2 . (3) 

Even though there is no actual physical system associated with the Gaussian model, it 
is commonly used in the literature as a toy model for nonlocal nonlinearities. Note that 
without loss of generality the width of the kernel K has been set to unity, in order to 
have the same scaling as used in [26"|l27]. 



3. Linear stability analysis of higher-order solitons 

Let $ be a bright solitonic solution to our governing equation flTJ 

d>M)=^ s (r)e iAt , (4) 

where A is the propagation constant or chemical potential for the case of optical beam 
or BEC, respectively, and denotes the stationary profile of the soliton. Because we 
will not consider solitons carrying angular momenta (e.g., vortices), we can choose ips{ r ) 
to be real. 

In order to find numerically exact stationary profiles ijjs{ r ), we use variational 
solutions as input to an iterative solver [30]. Typically, we use a grid of 400 x 400 points 
to determine ^s( r )- This transverse resolution is also employed for numerical integration 
of Eq. (pQ), i.e., for beam propagation or time evolution of the two-dimensional BEC. 

Figure [2] shows solitonic family curves or the two higher order solitons we choose 
to study here, the second-order radial state R 2 and the quadrupole Q. Apart from 
the total angular momentum, there are two conserved functionals, i.e. the Hamiltonian 
W[ip\ associated with invariance with respect to time-translations and the mass M[ip] 
due to a global U(l) phase-invariance: 

H\M = J |Ws| 2 d 2 r-i J |^ s (r)| 2 ir(r-r')|^(r')rd 2 r'd 2 r, (5) 
M[H = [ |^s| 2 d 2 r. (6) 



Obviously, the family curves for the R 2 and Q solitons are quite close to each other, 
which was used in [26] to explain the observed quasiperiodic shape transformations 
(energy crossing). However, we will see in the following analysis of projected propagation 
dynamics in Sec. HJthat this very intuitive picture does not hold. 

Let us first recall that the linear stability of solitonic solutions can be studied as 
an eigenvalue problem as follows. We introduce a small perturbation 5ifj(r,t) to our 
solitonic solution ips{ T ) vi & 

^(r,t) = [^ s (r) + ^(r,t)]e iAi (7) 

Plugging Eq. ([7J into the governing equation Eq. ([1]) and retaining only first order terms 
in Sip, yields the following (linear) evolution equation for 8ip: 



id t -\ + A + J K(\r - r'D^KOdr' 
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Figure 2. Solitonic family curves for the second-order radial soliton R2 (blue) and 
the quadrupole soliton Q (red). Dashed lines indicate parameter domains where the 
soliton is linearly unstable. 



+^s(r) J K(\v-v'\)^(v') [#(r',t) + ^(r',t)]dV = 0. 

With the ansatz 

5ip{r,t) = 5u{r)e iKt + 5v*(r)e- iKH 
for the perturbation we can derive the eigenvalue problem (BdG equation) 

A-A + J K(\r - r'DVKOdV 5u{r) 
+ip s (r) J K{\r - r'\)ip s {r') [Su{r') + 5v(r')] dV = k5u{y) 



(8) 
(9) 



A - A + 



j X(|r-r'|)Vi(r')dV 



Sv(t) 



-^s(r) J K(\r - r'\)^j s (r') [6v(r') + Su(r')) d 2 r' = kSv(t) 



(10) 



(11) 



Real-valued eigenvalues k of Eq. ( TTUj) are termed orbitally stable and the corresponding 
eigenvector (Su, Sv) can be chosen real-valued. On the other hand, complex eigenvalues 
with negative imaginary part indicate exponentially growing instabilities. We note that 
due to the special structure of Eq. ([TO]) [which has its origins in the Hamiltonian structure 
of Eq. ([[])], if k is an eigenvalue, then — n as well as ±k* are also eigenvalues. 

Next, we solve Eq. (|T0|) in order to obtain the internal modes of the second-order 
radial soliton R2 and the quadrupole Q, respectively. A trivial solution to this problem 
is always given by (5u,5v) = i(V'S) — ^s) with eigenvalue k — 0. This so-called trivial 
phase mode is linked to the phase invariance of solitons. Derivatives of this trivial phase 
mode with respect to x or y are also trivial eigenvector^] with eigenvalue k = 0, and thus 
the eigenvalue k = is degenerate. Moreover, due to symmetry properties of the system 
trivial modes appear twice in the spectrum, i.e., we expect sixfold degeneracy of the 
eigenvalue k = 0. However, when solving the discretized version of Eq. ( 1T0|) numerically, 
this degeneracy may be lifted. Thus, degenerate eigenvectors with zero eigenvalues may 



J The trivial modes (Su, 5v) — ±(d x ips,—d x i/)s) resp. (Su,6v) 
translational invariance of the system. 



±{dyips, —dyips) are linked to the 
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in fact become slightly complex without actually indicating an instability. In other 
words, their nonzero imaginary part is a numerical artefact of the discretization and 
occurs because the full eigenspace has to be spanned by the eigenvectors. The actual 
computation of the linear eigenvalue problem Eq. ( ITU]) is numerically expensive, since 
the matrix we have to diagonalize is full, i.e. all entries are nonzero. In order to achieve 
reasonable computation times, we usually reduce the grid-size to 100 x 100 points only. 
Then, the matrix we have to diagonalize has 4 x 10 8 non-zero elements. 

In Fig. [3l we show the spectrum of the linear stability analysis (BdG equation) for 
the second-order radial soliton R 2 and the quadrupole soliton Q [a) resp. b)] for mass 
M = 200. Note that for modes with purely imaginary eigenvalue k = ilmn, Eq. ([9]) 
reads 5ip(r,t) = [6u(r) + 6v*(r)] e ~ Im(K )' ; and it makes sense to define 



Only the second order radial state i? 2 is unstable, and we name the two unstable 
internal modes ei, e 2 . The unstable modes e 1; e 2 ought to be degenerate for symmetry 
reasons, the small splitting of the eigenvalues («i ~ —2.7i, k 2 « — 2.5i) is again 
a numerical artefact due to the discretization of the eigenvalue problem Eq. ([TO]) . 
Interestingly, the shape of the unstable eigenmodes ei(r), e 2 (r) resembles quadrupoles. 
In fact, for practical purposes (see next section) as well as to verify these findings we 
furthermore solved the eigenvalue problem Eq. (jTOl) for R 2 on a radial grid [31] with 
eightfold resolution. Then, instead of two stable and unstable quadrupoles, one finds 
one stable and unstable vortex with topological charge m = ±2 and |«| ~ 2.74. The 
vortices corresponding to m = 2 and m = — 2 can be superposed to again yield the 
quadrupoles e±, e 2 found already with the full 2D solver, but with much higher precision. 
Because Eq. ( ITU]) is linear, the amplitudes of the ij are not fixed, and we normalize the 
latter according to 



The quadrupole soliton Q in Fig. [3b) is stable, because all complex eigenvalues 
correspond to trivial modes and hence the complex form of these eigenvalues is a 
numerical artefact as discussed above. However, the quadrupole becomes linearly 
unstable for M < 90, as indicated in Fig. [2] by dashed lines. In Fig. HJ we show 
the results of our numerical stablility analysis for the quadrupole soliton Q with mass 
M = 85. Interestingly, the unstable mode e.\ with K\ ~ —1.2i resembles the second- 
order radial soliton i? 2 , i.e., a hump with a (modulated, i.e. not rotationally symmetric) 
ring. 

4. Projected nonlinear dynamics 

The typical dynamics for i? 2 (here for M = 200) as an initial condition is shown in 
Fig. [5] a). To determine the shape of i? 2 , we use the iterative solver mentioned above on 
a grid containing 400 x 400 points, and we use the same grid for the actual propagation. 



e(r) = 8u(r) + 5v*(r). 



(12) 




(13) 



Quasiperiodic oscillations and homoclinic orbits in the NNLS 



a) 



trivial 
/modes 



a 



-3L 




-'20 



o " 



20 -20 



Q 


b) 




trivial 




y modes 




■f ■■ 







20 



Figure 3. Spectrum of the linear stability analysis (BdG equation) centered around 
zero for a) the second-order radial soliton i?2, and b) the quadrupole Q. Both solitons 
have mass M = 200. The radial soliton Ri exhibits instabilities and the unstable 
eigenmodes ei(r), (=2 (r) resemble quadrupoles [see two insets in a)]; the quadrupole 
soliton Q is stable. For both solitons, the degeneracy of the trivial modes is lifted, a 
numerical artefact due to the discretization of the eigenvalue problem Eq. ([TO]) . For 
sake of readability, the insets in a) show the absolute square |e(r)| 2 = |<5it(r) +<5w*(r)| 2 
only. 
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Figure 4. Spectrum of the linear stability analysis (BdG equation) centered around 
zero for the quadrupole Q with mass M = 85. The unstable eigenmode ei(r) resembles 
the shape of R2 (see inset), but is of course not rotationally symmetric. Again, the 
degeneracy of the trivial modes is lifted, a numerical artefact due to the discretization 
of the eigenvalue problem Eq. (|10[) . For sake of readability, the inset shows the absolute 
square | e(r) | 2 = \5u(r) + Sv*(r)\ 2 only. 



As we have seen in Sec. [31 the second-order radial soliton R 2 is unstable over the whole 
range of mass M and therefore any perturbation, that has a non-zero overlap with 
the unstable internal modes e\, £2 will lead to an exponential growth of the latter. 
Practically, the residual in numerical determination of R2 as well as the propagation 
algorithm based on the Fourier split-step method [Tj lead to inevitable numerical noise 
when propagating and therefore trigger the instability without adding any additional 
perturbation. In our case, however, we added the eigenmode t\ as initial perturbation 
with tiny amplitude ~ 10 -4 to the soliton R2 to control the breakup in a preferred 
direction. For small times the dynamics is governed by the exponential growth of the 
unstable internal mode e\, while for later times the evolution becomes highly non-linear, 
exhibiting oscillations between R 2 [see inset (a) in Fig. [5] a)] and a state that resembles 
the quadrupole soliton Q [see inset (ft) in Fig. [5] a)] [26]. This state ((3) we will call the 




Figure 5. a) Evolution of the peak-intensity of the second-order radial soliton i?2 with 
mass M = 200 (upper blue curve). As expected from the stability analysis Fig. E]b), 
the peak-intensity of the quadrupole soliton Q with same mass (lower red curve) is 
constant during propagation. Figure b) shows the projected dynamics in the variables 
U(t),S(t) [see Eq. (fT8|)] for initial conditions R2 [blue curve, starting at (a)} and Q 
[red curve, starting at (7)]. The shape of the former curve hints at a homoclinic 
connection, where the homoclinic point corresponds to R2 {a). Figure c) presents 
the same dynamics as b), with an additional dimension given by the variable w [see 
Eq. (1191) ]. In this three-dimensional projection, the distance between the quadrupole 
Q (7) and the "turning point" (/?) becomes apparent. For reasons of clarity, the 3D- 
dynamics (blue) is additionally projected into (S, w)-plane (black), and the orbit of the 
the quadrupole is again shown in red. The three insets show snapshots of the nonlinear 
dynamics. 



"turning point". In the following we will examine in detail the origin and properties of 
these oscillations. 

4-1. Projection methods 

Let us now introduce the projection method mentioned in the introduction J28j|29] and 
adopt it to our problem. To this end, we recall the scalar product of two complex 
functions / and g, defined as 

(f,g) = f n*)g(r)dh. (14) 

Obviously, the (unstable) internal modes e 3 - of R2 introduced before [see Fig. [3] are not 
orthogonal to their complex conjugate (stable) e* (j = 1, 2) counterparts with respect 
to this inner product. In other words, stable and unstable eigenspaces E s and E u 
spanned by eigenf unctions {e*, e^} resp. {ei, e 2 } are not mutually orthogonal. Thus, 
straightforward projections onto e.j and e* do not help to elucidate the propagation 
dynamics. To overcome this difficulty we introduce a set of functions which is 
biorthogonal to ej,e* using a Gram-Schmidt-like technique as follows. First, we define 

e j± = ej - (e*,ej)e*, (15) 

which is simply the projection of the unstable eigenmode ej onto the orthogonal 
complement of the stable eigenmode e*. Second, we note that (ej±)* = e* — (ij,e*)ij 
corresponds to projection of the stable eigenmode e* onto the orthogonal complement 
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Figure 6. Schematic sketch of the relation between ej,e*,ej 



and {e j± )*. By 
is orthogonal 



construction, ej± is orthogonal to the stable eigenvector e*, and (ej±) 
to the unstable eigenvector Sj. It is worth to notice that ej± and (e^jj* are not 
orthogonal to each other, b) and c) show the modulus squared of the internal mode 
ei and eu_, respectively. 



of the unstable eigenmode e,-. Then, it is easy to verify biorthogonality of e,j_, (e 



with respect to e,j, e*: 



(ej,(e j± y) = 0(ej,e j± ) ^ 



J-LJ 



j±, y^j±) 
(16) 

is depicted, 



In Fig. [6^) a schematic sketch of the relation between ij, e*, ej±, and (e 
and Fig. [6b-c) show e.\ and e q ± explicitly. It is worth to notice that (ej±, (ej±)*) ^ 0, 
i.e., ej± and (ejj_)* are not orthogonal to each other. 

In order to analyze the propagation dynamics of a solution if)(x,t) of Eq. (JTJ), we 
introduce the quantities 

Uj = (e j± ,i>), S j =((e j± )*,i>). (17) 

By construction, Uj is associated with the unstable eigenmode only (ej± is orthogonal 
to the stable one), while Sj is associated with the stable eigenmode only. Finally, for 
i?2, the two unstable eigenvectors i\, e~2 are degenerate (due to rotational symmetry 
about the origin), therefore we introduce the rotationally invariant projected variables 



U(t) 



\ 



5>, 



S(t) 



\ 



Then, any pair of wavefunctions ipi(x, t) and ipzfe, t) related through a rotation amounts 



to the same value of U(t) and S(t). For a rigorous proof, see Appendix Appendix A 



4-2. Indication of homoclinic connections 

Figure [5] b) illustrates the dynamics shown in Fig. |5] a) in the variables S(t),U(t) 
introduced in Eq. ( fl8l) . We clearly see the second-order radial soliton R 2 (a) decaying 
into a quadrupole-like state (/3), the "turning point", and then coming back to R 2 - In 
the vicinity of R 2 , the decay starts via the local unstable eigenspace E u (i.e., U(t) > 0, 
S(t) w 0), and the revival of R 2 happens via the local stable eigenspace E s (i.e., S(t) > 0, 
U(t) ~ 0). The fact that the system repeatedly returns (close) to its initial state R 2 and 
remains at this point some finite, non-constant time with (nearly) zero velocity, hints at 
the existence of a homoclinic connection. A homoclinic connection is a solution which 
is asymptotic to R 2 both in the t — > oo and t — > — oo limit. The time-span, in which the 
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solution remains close to its initial state R 2 , i.e. the homoclinic point (a) in Fig. [5]b), 
with practically zero velocity, corresponds to intervals with maximum (nearly) constant 
peak-intensities in Fig. |5]a). Because we added a small perturbation in the direction 
of the eigenmode e.\ to the initial condition R 2 , and the presence of numerical noise in 
general, we do not see the exact homoclinic connection in our numerical simulations; as 
the trajectory comes back towards R 2 along E s , there is always a small perturbation 
along the unstable eigenspace E u and the trajectory leaves the neighborhood of R2 
to return to it later on. We want to stress here that the existence of homoclinic 
connections is by no means anticipated in general; our numerical results however indicate 
the existence of such homoclinic connections and their persistence along a large range 
of the mass M. 

To further illustrate that the "turning point" (f3) is indeed well-separated from the 
quadrupole soliton Q (7), we introduce a third variable w by projecting the solitonic 
wave function ip onto the radial soliton R 2 , 

Obviously, for ip = R 2 we find w = 1, while for if) = Q for symmetry reasons we have 
w — 0. Figure [5] c) shows the resulting projected dynamics on the variables U,S,w. 
We clearly recognize similarities with Fig. |5]b), however, it becomes much more clear 
how the solution evolves from its origin (a) and becomes much more "quadrupole-like" 
in {13). In particular, the important separation between the quadrupole-like "turning- 
point" (/?), which still maintains a nonzero projection on R 2 and the quadrupole soliton 
Q (7) becomes evident. 



4-3. Quasiperiodic motion 

In the previous section, we argued that due to numerical limitations, we cannot actually 
track the homoclinic orbit precisely, but what we find are trajectories that are very close 
to the homoclinic connection. In the present section we will further probe the dynamical 
importance of the homoclinic orbit by studying trajectories adjacent to it. In a sense, 
the "turning point" (/3) of the homoclinic orbit is a state "in between" R 2 and Q. Here 
we will investigate the dynamics of such "in between" states obtained by perturbing the 
homoclinic orbit at the "turning point" ((3). The perturbations we will consider are not 
necessarily small and, as we will see, they typically lead to quasiperiodic oscillations. 

A homoclinic orbit is obtained by (slightly) perturbing the initial wavefunction of R 2 
in the direction of one of the unstable modes (e.g. of ei) and integrating Eq. ([1]) forward 
in time. Choosing the direction of the initial perturbation fixes the "orientation" of the 
subsequent dynamics, and we can thus decompose the wavefunction at the turning point 
[point ((3) in Fig. |5] t t into a part parallel to the quadrupole soliton Q and a remainder 
L 



4>(r,t t ) = c Q Q(r)+L(r), 



(20) 
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where cq = (Q, ip(r,t t ))/{Q, Q) was introduceq§|. Perturbed wavefunctions Vr( r ) are 
then constructed through 

^(r) = c Q Q(r) + r L(r) (21) 

where V parametrizes mixed states between R2 and Q, and, in Eq. (1221) . the wavefunction 
was normalized. Clearly, for T = 1, the homoclinic trajectory of R2 can be recovered, 
whereas of F = 0, the quadrupole soliton is recovered. In the following the time evolution 
of the function ipr will be studied. 

Let us first consider the dynamics for T = 1.01 as shown in Fig. [71 which indicates 
quasiperiodic behavior for small times (up to t ~ 25). The time spend by this orbit close 
to i?2 is much smaller than for the homoclinic connection, of the previous section. This 
becomes apparent when comparing the peak-intensity evolution in Fig. [7h) with the one 
in Fig. [5^). In the (U, S, w) projection, this fact results in a smoother curve close to the 
origin (whereas for a homoclinic connection a kink appears as R2 is approached, while 
the "velocity" approaches zero). On the other hand, in the intensity representation, 
Fig. H>h), the difference between homoclinic and quasiperiodic behavior is much harder 
to discern. Propagation in Fig. [7Ji) and b) is shown until t = 35, when the dynamics 
already deviates from the quasiperiodic orbit, indicating that the latter is unstable. This 
behavior hints to the existence of some chaotic region in state-space, an issue that will 
be studied elsewhere. 

On the other hand, the dynamics for T = 0.99, shown in Fig. [HI appear again 
quasiperiodic (see also the discussion of the Fourier spectra in Sec. I4.4p . but in this case 
the orbit appears stable, as it persists at least up to t = 1500. The qualitatively different 
behaviour for T = 1.01 and T = 0.99 with respect to stability further corroborates the 
importance of the homoclinic solution T = 1.00 (i^)- in a certain sense, the homoclinic 
orbit "organizes" regions of stability in parameter space. However, the homoclinic 
orbit should not be seen as a kind of "boundary" between regions of different stability 
behaviours, because it is just a one-dimensional line in the highly- dimensional parameter 
space. 

Let us finally consider the trajectory in Fig. [9] which is far away from both the 
quadrupole soliton as well as from the "turning point" (/3) by letting T = 0.5. The 
dynamics is still quasiperiodic and stable (at least up to t = 1500), but involves multiple 
frequencies. Interestingly, the dominant frequency of oscillation with period T « 2.6 
can be related to a stable eigenvalue of the quadrupole soliton Q for M = 200. In the 
(stable) eigenvalue spectrum of Q shown in Fig. [3b), the internal mode with n rs 2.6 
resembles a (modulated) ring with a hump (not shown). The duration of one period 
T would then be given by T = 2tt/k 2.4, which is what we find when we slightly 

§ More generally, if the direction of the breakup is arbitrary, one may generalize Eq. (|20jl by 
decomposing ^(it) into two quadrupoles Q\, Q27 where Q\ is rotated by it /A with respect to Q2, 
via V(r,t t ) = c Ql Qi(r) + c Q2 Q 2 (r) + L(r). 
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Figure 7. Evolution of ipr for T = 1.01 denned in Eq. (f2~Tj) . (a) shows the peak- 
intensity, (b) the orbit in lower-dimensional S, U, w representation, and (c-h) snapshots 
of the dynamics. The coloring is the same as in Fig. [5l where the blue curve again 
represents the actual 3D dynamics and the black curve its projection on the (s,w)- 
plane, and the red curve is the orbit of the quadrupole. 
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Figure 8. Evolution of for T = 0.99 defined in Eq. ([2~T|) . (a) shows the peak- 
intensity, (b) the orbit in lower-dimensional S, U, w representation and (c-h) snapshots 
of the dynamics. The coloring is the same as in Fig. [5l where the blue curve again 
represents the actual 3D dynamics and the black curve its projection on the (s,w)- 
plane, and the red curve is the orbit of the quadrupole. 



perturb the quadrupole soliton Q by this mode. Moreover, for T = 0.1 (not shown) we 
also find an oscillation with period T rs 2.4. In both case, the propagation dynamics 
resemble the one shown in Fig. |9]for T = 0.5. Thus, even though for T = 0.5 we are 
no longer in the region where perturbation analysis of the quadrupole soliton Q holds, 
we still find qualitatively similar dynamics. We note that in the same system Eq. ([T|), 
quasiperiodic nonlinear solutions (so-called azimuthons) linked to stable internal modes 
of solitons were reported earlier [3lJ[32] . 

To sum up, we have identified a family of stable quasiperiodic solutions to Eq. ([T|), 
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Figure 9. Evolution of for T = 0.5 denned in Eq. (|2Tj) . (a) shows the peak- 
intensity, (b) the orbit in lower-dimensional S, U, w representation and (c-h) snapshots 
of the dynamics. The coloring is the same as in Fig. [5l where the blue curve again 
represents the actual 3D dynamics and the black curve its projection on the (s, w)- 
plane, and the red curve is the orbit of the quadrupole. 



starting from ijj-p given in Eq. (12 ip and < V < 1. The two limiting solutions are the 
stable quadrupole solitons Q (T = 0) and the homoclinic orbit linked to the unstable 
radial solitons R2 (T = 1). We want to emphasize here that for lower masses, where the 
quadrupole soliton Q becomes unstable (e.g., M = 85), we were not able to find stable 
quasiperiodic solutions by the same construction. 



4-4- Fourier spectrum 

Further insight can be gained by considering the Fourier spectrum of the above 
trajectories. Given a trajectory ifj(r, t) we compute the modulus of the Fourier transform 
J 7 of the wavefunction at a fixed point in space (in our case the origin r = 0): 

/(a;) = |^(r = 0,t)| 2 . (23) 

For a bright soliton solution of the form Eq. one would expect f(oo) to comprise of 
a single sharp peak at u = A. On the other hand, in the case of quasiperiodic dynamics 
with vibration frequency Q and propagation constant A, one would expect peaks at 
A + iriQ, where m is integer. This is readily verified for the orbits with a = 0.99 and 
a = 0.5, as can be seen in Fig. UHl where we see sharp peaks associated with these orbits. 
On the other hand, there is no well defined periodicity associated with the homoclinic 
orbit, since the time spent in the vicinity of R2 is in principle infinite. In practice, 
this time is greatly affected by numerical noise and the spectrum appears continuous 
[see Fig. HDk)]. Even if it is possible to associate a dominant frequency Q with the 
homoclinic orbit, f{oo) around Q is much broader than in the case of quasiperiodic 
orbits for T = 0.99 and T = 0.5 [see Fig. [10b) and c)jjj 

|| A limitation on the spectral resolution for /(w) for the homoclinic orbit appears due to the fact 
that dynamics become unstable around t — 520. Here, we used the interval t = [0 : 500] to compute 



Quasiperiodic oscillations and homoclinic orbits in the NNLS 



14 



Thus, the Fourier spectra yield an additional indication of the qualitatively different 
nature of the dynamics of Sec. 14.21 from the quasiperiodic motion of Sec. I4.3[ providing 
further support for the conjectured existence of an underlying homoclinic connection in 
the former case. 

5. Conclusions 

In previous works, an oscillatory shape-transformation of modes in nonlocal media has 
been observed [261 [27]. m this paper, we approached this phenomenon by means of 
linear stability analysis and projection techniques borrowed from dynamical systems 
studies of dissipative PDEs. By studying the linear stability of the quadrupole soliton 
Q and the second-order radial soliton R 2 , we found that the former becomes linearly 
stable for mass M > 90, whereas the latter remains linearly unstable for all masses. 
The initial stage of the shape-transformations under consideration, i.e. the emergence 
of a new state on top of R 2 , can be understood in terms of this linear instability, which 
is triggered by the unavoidable numerical noise. However, the most striking feature of 
the dynamics, i.e. the return to the initial state, is inherently nonlinear, as it occurs 
only after the linear instability saturates. To study this phenomenon, we introduced a 
low-dimensional representation of the dynamics, through a projection to dynamically 
important states, which were constructed from the radial soliton R 2 itself and its 
unstable/stable eigenmodes. Projecting the time evolution of the wavefunction ip(r,t) 
(obtained by integrating the NLS) onto these states allows a visualization of oscillatory 
shape-transformations in terms of trajectories, revealing that shape-transformations 
can be interpreted as a homoclinic orbit leaving and re-approaching R 2 . Moreover, 
in the neighborhood of this homoclinic orbit we found quasiperiodic solutions, which 
for small enough perturbations resemble the homoclinic connection. This indicates 
that the homoclinic connection provides a basic recurrence mechanism around which 
quasiperiodic dynamics is organized, as is common in lower-dimensional dynamical 
systems [33]. We were also able to construct and identify a whole family of stable 
quasiperiodic orbits when the quadrupole soliton Q is stable. 

The projection method introduced here allows a compact representation of the 
dynamics, dual to the commonly used intensity plots. Moreover, in certain cases it 
helps to uncover features of the dynamics that are not apparent in snapshots of the 
intensity evolution. We expect that similar studies can be carried out for other states 
exhibiting similar dynamics [26] and that our projection method (or similar extensions 
of the methods of Refs. [28],|29]) could be applied to a variety of high- and infinite- 
dimensional conservative systems. 

the spectrum. Thus, compared to the other two spectra shown in Fig. 1 1 01 where the propagation was 
performed until t — 1500, the spectral resolution is coarser by a factor of three. 
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Figure 10. a) Spectrum f(to) — \F{4'{ Y — 0, i ) | 2 corresponding to the homoclinic 
orbit r = 1.00 (red), and quasiperiodic orbits with T — 0.99 (black) and T — 0.5 (blue) 
in logarithmic scale, b) Same information in linear scale. c)-e) show magnifications of 
single peaks of b). 
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Appendix 

Appendix A. Rotational invariance of U(t) and S(t) 

Here we prove that the quantities U(t), S(t) are rotationally invariant, i.e. they have the 
same value if we substitute i[)(x,y,t) with 1Z(9)i[)(x,y,t) = ip(x cos 9 — ysin9,x sin 9 + 
ycos9,t), where 11(9) is an SO(2) rotation. 

The eigenproblem Eq. ( ITUl) for the ring soliton R2 is rotationally symmetric and, 
as a result, its internal modes e\, e 2 transform according to 

2 

H(6)e i = ^2D ji {e)e j , (A.l) 

j=i 

where D(#) is a two-dimensional matrix-representation of SO(2). The explicit 
representation D(9) depends on the basis ij, but for our purposes it is sufficient to 
show that we have a real representation. We begin by noting that the constraints of 
orthogonality, D T D = 1, and unit determinant, det(D) = 1, lead to the following general 
form 

DM - ( a{6) m \ (A ?) 

{ ) ~ { -P*(0) a*(9) ) (A ' 2) 

where the functions a(9), (3(9) are related through 

det (D(0)) = \a(9)\ 2 + \(3(9)\ 2 = 1. (A.3) 

On the other hand, using e 2 = H(9 )ii, where 9 is the angle that rotates e± onto e 2 , 
we can express all matrix elements Djj = (ej,1Z(9)ei) in terms of Dn, 

mm I a ^ ^(9 + 9 ) , 

m ~ 1 a(9- 9 ) a(9) ' (A ' 4) 
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Comparing with Eq. (I A . 2 [) we conclude that a(6) = a* (8) and thus our representation 
is real, and that a(6 — 6 ) = —a(6 + ^o)0 

Using Eqs. ( lA.ll) -( fA~3l) in definition Eq. ( 1T5|) . along with the relation (e*,ei) = 
(e^, e 2 ), one can show that 

2 

n{e)e i± = ^2D ji {6)e j± . (A.5) 

Then, using Eqs. ( 1A.2l) -( )A~5l) . it's easy to show that 

U 2 (t) = | (e 1± ,K(6)iP) | 2 + | (e 2± , ft W> | 2 

= |(^M)e 1± ,^)| 2 + |(^M)e 2± ,^)| 2 
= |(eix,^)| 2 + |(e 2± ,^)| 2 

= ^) ■ 
A similar proof holds for S(t). 



% In our numerical results 9q — 7r /4 and one can see that our representation is in fact equivalent to 

/ cos(2tf) -sin(20) \ 
1 j V sin(20) cos(2#) y ' 



